Os dados

set.seed(12345)

lastfm = read_csv(here::here("data/experimento-lastfm.csv"), 
                  col_types = cols(.default = col_double(), 
                                   user = col_character()))

lastfm = lastfm %>% 
  sample_n(300) %>% 
  select(news, old, mediana_pop)

glimpse(lastfm)
## Observations: 300
## Variables: 3
## $ news        <dbl> 28, 35, 13, 24, 14, 17, 13, 21, 34, 55, 10, 33, 10, …
## $ old         <dbl> 61, 194, 70, 96, 130, 67, 106, 123, 76, 78, 76, 116,…
## $ mediana_pop <dbl> 6.105585, 5.376812, 5.713082, 4.564335, 5.782320, 5.…

Proporção de artistas novos e popularidade

Utilizaremos ICs para estimar duas métricas sobre os usuários do LastFM em geral durante um período de 6 meses. Em ambos os casos faremos isso a partir de uma amostra de 300 usuários. As duas métricas são: 1. Qual a proporção de novos artistas em geral escutada por usuários? 2. Para os usuários que gostam de música muito pop (mediana_pop > 5), qual a correlação entre a popularidade mediana dos artistas escutado e a proporção dos artistas escutados que eram novos.

Perguntas

Os questionamentos levantados serão respondidos nos itens abaixo

1. Qual a proporção de novos artistas em geral escutada por usuários?

A média da proporção de novos artistas em geral escutada pelos usuários é dada pela variável theta_c:

set.seed(9098)

theta_c = lastfm %>% 
    mutate(news_prop = news/(news + old)) %>%
    pull(news_prop) %>%
    mean()

theta_c
## [1] 0.2483568

Aplicando a técnica de bootstrap (com 3000 repetições) para calcular a média da proporção de artistas novos de amostras aleatórias retiradas da amostra principal, temos:

set.seed(9098)

repeticoes = 3000

bootstrap_q1 <- function(ds) {
  ds = ds %>% mutate(news_prop = news/(news + old))
  news_prop = ds %>% pull(news_prop)
  bootstrap <- sample(news_prop, size = NROW(news_prop), replace = TRUE)
  
  return(mean(bootstrap))
}

reamostragens = tibble(i = 1:repeticoes) %>% 
  mutate(theta_c_s = map_dbl(i, ~ bootstrap_q1(lastfm)))

Agora, vamos visualizar as reamostragens que obtivemos:

set.seed(9098)

graf_reamostragens_q1 = reamostragens %>%
  ggplot(aes(x = theta_c_s)) +
  geom_histogram(binwidth = .00125,
                 fill = "pink",
                 colour = "black")

ggplotly(graf_reamostragens_q1)

Pela visualização acima, além de percebermos que o gráfico segue uma distribuição normal, vemos que a maioria das médias das reamostras estão entre 24 e 26%.

Agora, calculando o intervalo de confiança (IC):

set.seed(9098)

ic_q1 = reamostragens %>% 
  mutate(erro = theta_c_s - theta_c) %>% 
  summarise(erro_i = quantile(erro, .05), 
            erro_s = quantile(erro, .95))

ic_q1 = ic_q1 %>% 
  mutate(valor_i = theta_c + erro_i, 
         valor_s = theta_c + erro_s)


ic_q1

Agora, vamos visualizar o IC (obtido com a técnica de bootstrap):

set.seed(9098)

graf_ic_q1 = ggplot() +
  geom_rect(
    data = ic_q1,
    aes(xmin = valor_i, xmax = valor_s),
    ymin = -Inf,
    ymax = Inf,
    fill = "lightgreen",
    alpha = .25
  ) +
  geom_histogram(
    data = reamostragens,
    aes(theta_c_s),
    binwidth = .0015,
    fill = "pink",
    colour = "black"
  ) +
  geom_vline(xintercept = theta_c, color = "blue")

graf_ic_q1

Dessa forma, podemos afirmar com 95% de confiança que, com uma amostra de 300 itens, a média populacional da proporção de novos artistas escutados pelos usuários será estimada entre 23,76 e 25,97%.

Agora, vamos comparar o resultado obtido com um bootstrapper já implementado pela biblioteca boot:

set.seed(9098)

theta <- function(d, i) {
    d = d %>% slice(i) %>%
        mutate(news_prop = news/(news + old)) %>% 
        summarise(media = mean(news_prop))
    
    m = d %>% pull(media)
    m
}

booted = boot(data = lastfm,
              statistic = theta,
              R = 3000)

ci = tidy(booted, 
          conf.level = .95,
          conf.method = "bca",
          conf.int = TRUE)

glimpse(ci)
## Observations: 1
## Variables: 5
## $ statistic <dbl> 0.2483568
## $ bias      <dbl> 0.0001609142
## $ std.error <dbl> 0.006733435
## $ conf.low  <dbl> 0.2356069
## $ conf.high <dbl> 0.2616769

Utilizando um bootstrapper já implementado podemos ver que, com 95% de confiança, a média populacional da proporção de novos artistas escutados pelos usuários será estimada entre 23,56 e 26,16%. Esse intervalo obtido com um bootstrapper já implementado condiz com o que obtivemos previamente.

2. Para os usuários que gostam de música muito pop (mediana_pop > 5), qual a correlação entre a popularidade mediana dos artistas escutado e a proporção dos artistas escutados que eram novos.

Primeiramente, filtrando nossos dados (apenas usuários com mediana_pop > 5):

lastfm_q2 = lastfm %>% filter(mediana_pop > 5)
theta_q2 <- function(ds) {
    ds = ds %>% mutate(news_prop = news/(news + old))
    
    c <- cor(ds$news_prop, ds$mediana_pop)
    return(c)
}    

theta_c = theta_q2(lastfm_q2)

theta_c
## [1] -0.088961

Como o valor obtido é próximo de zero, aparentemente não há nenhuma correlação.

Vamos agora visualizar um gráfico correlacionando as variáveis:

lastfm_q2 %>% mutate(news_prop = news/(news + old)) %>%
    ggplot(aes(x = news_prop, y = mediana_pop)) + 
    geom_point()

Vemos pontos muito dispersos, a correlação entre a popularidade mediana dos artistas escutado e a proporção dos artistas escutados que eram novos é baixa.

Realizando o bootstrapping:

set.seed(9098)

repeticoes = 3000

bootstrap_q2 <- function(ds) {
  ds = ds %>% mutate(news_prop = news/(news + old))
  
  bootstrap <- sample_n(ds, 
                        size = NROW(news_prop), 
                        replace = TRUE)
  
  return(cor(bootstrap$news_prop, bootstrap$mediana_pop))
}

reamostragens = tibble(i = 1:repeticoes) %>% 
    mutate(theta_c_s = map_dbl(i, ~ bootstrap_q2(lastfm_q2)))

graf_reamostragens_q2 = reamostragens %>%
  ggplot(aes(x = theta_c_s)) +
  geom_histogram(binwidth = .01,
                 colour = "black",
                 fill = "pink")

ggplotly(graf_reamostragens_q2)

Como há menos elementos em cada classe ocorreram mais resultados diferentes quando comparamos com a visualização das reamostragens feitas na questão anterior. Percebe-se também que a correlação é próxima bem próxima de 0 em todas as reamostragens.

Calculando o IC:

set.seed(9098)

ic_q2 = reamostragens %>% 
  mutate(erro = theta_c_s - theta_c) %>% 
  summarise(erro_i = quantile(erro, .05), 
            erro_s = quantile(erro, .95))

ic_q2 = ic_q2 %>% 
  mutate(valor_i = theta_c + erro_i, 
         valor_s = theta_c + erro_s)


ic_q2

Visualizando esse IC:

graf_ic_q2 = ggplot() +
  geom_rect(
    data = ic_q2,
    aes(xmin = valor_i, xmax = valor_s),
    ymin = -Inf,
    ymax = Inf,
    fill = "lightgreen",
    alpha = .25
  ) +
  geom_histogram(
    data = reamostragens,
    aes(theta_c_s),
    binwidth = .01,
    fill = "pink",
    colour = "black"
  ) +
  geom_vline(xintercept = theta_c, color = "blue")

graf_ic_q2

Observamos com 95% de confiança que, a correlação entre a popularidade mediana dos artistas escutado e a proporção dos artistas escutados que eram novos está estimada entre -0,2 e 0,026. Como esse intervalo contém o valor 0, pode-se afirmar que, caso a correlação exista (pode acontecer da correlação ser zero), ela é bem próximo de zero.

Agora, vamos comparar utilizando o bootstrapper da biblioteca boot:

set.seed(9098)

theta <- function(d, i) {
    d = d %>% slice(i) %>%
        mutate(news_prop = news/(news + old)) %>% 
        summarise(cor = cor(news_prop, mediana_pop))
    
    c = d %>% pull(cor)
    c
}

booted = boot(data = lastfm_q2,
              statistic = theta,
              R = 3000)

ci = tidy(booted, 
          conf.level = .95,
          conf.method = "bca",
          conf.int = TRUE)

glimpse(ci)
## Observations: 1
## Variables: 5
## $ statistic <dbl> -0.088961
## $ bias      <dbl> 0.002479051
## $ std.error <dbl> 0.06877761
## $ conf.low  <dbl> -0.2259757
## $ conf.high <dbl> 0.04092692

Com um bootstrapper já implementado, obtemos que, com 95% de confiança, a correlação entre a popularidade mediana dos artistas escutado e a proporção dos artistas escutados que eram novos está estimada entre -0,23 e 0,041%. Novamente percebemos que o zero está contido no intervalo, o que indica que a correlação pode não existir ou ser muito baixa.